Loading Data

Set wd and load packages.

As stars proxy

Intensity data is handled here.

Spatial Extent

The desired spatial extent can be set via subsetting: [extent_object].

Time Dimension

dates <- c(ymd("2018-01-03"), ymd("2018-01-15"), ymd("2018-01-27"), ymd("2018-02-08"), ymd("2018-02-20"), ymd("2018-03-04"), ymd("2018-03-16"), ymd("2018-03-28"), ymd("2018-04-09"), ymd("2018-04-21"), ymd("2018-05-03"), ymd("2018-05-15"), ymd("2018-05-27"), ymd("2018-06-08"), ymd("2018-06-20"), ymd("2018-07-02"), ymd("2018-07-14"), ymd("2018-07-26"), ymd("2018-08-07"), ymd("2018-08-19"), ymd("2018-08-31"), ymd("2018-09-12"), ymd("2018-09-24"), ymd("2018-10-06"), ymd("2018-10-18"), ymd("2018-10-30"), ymd("2018-11-11"), ymd("2018-11-23"), ymd("2018-12-05"), ymd("2018-12-17"), ymd("2018-12-29"))
int_vv <- st_set_dimensions(int_vv, 3, values = dates, names = "time")
int_vh <- st_set_dimensions(int_vh, 3, values = dates, names = "time")

Rain

Needs to be freshly warped for some reason, can’t be saved (write_stars) and then loaded. But: rain_warped.tif is a file with already aggregated rain data (31 time steps instead of 365).

Combine

Then combine all three. Give along = "bands" so that all attributes are treated as bands.

switch bands to attributes

This is more of a design choice.

comb_split <- split(comb, "bands")
names(comb_split) <- c("VV", "VH", "prec")
comb_split
## stars object with 3 dimensions and 3 attributes
## attribute(s):
##       VV               VH              prec        
##  Min.   :-20.54   Min.   :-24.74   Min.   : 0.00   
##  1st Qu.:-15.68   1st Qu.:-20.72   1st Qu.: 0.00   
##  Median :-14.50   Median :-19.72   Median : 0.20   
##  Mean   :-14.17   Mean   :-19.59   Mean   : 2.59   
##  3rd Qu.:-12.66   3rd Qu.:-18.45   3rd Qu.: 3.40   
##  Max.   : -5.13   Max.   :-10.62   Max.   :17.60   
##  NA's   :162998   NA's   :162998   NA's   :162998  
## dimension(s):
##      from   to     offset   delta                refsys point values x/y
## x    4745 4865     715082      10 ETRS89 / UTM zone 32N FALSE   NULL [x]
## y    1083 1202    5864854     -10 ETRS89 / UTM zone 32N FALSE   NULL [y]
## time    1   31 2018-01-03 12 days                  Date    NA   NULL

Data Exploration

Time Series for One Single Pixel

# get one pixel ts (at 10/10)
pixel.df <- as.data.frame(comb_split[1,100,100,])
pixel.df <- cbind(pixel.df, comb_split[2,100,100])
pixel.df <- cbind(pixel.df, comb_split[3,100,100])
pixel.df <- pixel.df[,c(3,4,8,12)]
# calculate additional measures
names(pixel.df) <- c("time", "VV", "VH", "prec")
pixel.df$sum <- pixel.df$VV + pixel.df$VH
pixel.df$diff <- pixel.df$VV - pixel.df$VH
pixel.df$ratio <- (pixel.df$VV / pixel.df$VH)
# plot all
ggplot() + 
  geom_bar(aes(x=pixel.df[,1], y=pixel.df$prec), stat='identity') +
  geom_line(aes(x=pixel.df[,1], y=pixel.df$VV, color="VV")) + 
  geom_line(aes(x=pixel.df[,1], y=pixel.df$VH, color="VH")) + 
  geom_line(aes(x=pixel.df[,1], y=pixel.df$sum, color="sum")) + 
  geom_line(aes(x=pixel.df[,1], y=pixel.df$diff, color="diff")) + 
  geom_line(aes(x=pixel.df[,1], y=pixel.df$ratio, color="ratio")) +
  scale_color_manual(name="Polarization", values=c(VV="red", VH="blue", sum="lightgreen", diff="forestgreen", ratio="yellow"))

# take a look at the ratio
ggplot() +
    geom_line(aes(x=pixel.df[,1], y=pixel.df$ratio))

A decrease in VV together with not so much decrease in VH might point to flooding with intact vegetation. More information on the plant appearance is needed.

Time Series for all (a Number of) Pixels

Giving the opportunity of constructing very different results. Depending on whether neighboring pixels or random pixels are selected, similarities and differences are seen.

# make first attribute, pixel 10/10 into dataframe template
pix.df <- as.data.frame(comb_split[1,100,100,])[,3:4]
# make data frame for all pixels
for (i in 1:120) {
  for (j in 1:119) {
    if (is.na(comb_split[1,i,j,1][[1]])) {
    }
    else {
      pix.df <- cbind(pix.df, as.vector(comb_split[1,i,j,][[1]]))
    }
  }
}
# delete template data frame column
pix.df <- pix.df[,c(1,3:9264)]
# assign numbers as names
names(pix.df) <- c("time", seq(1,9262,1))

Here we have the possibility to select which pixels are plotted. For example each 37th pixel in the dataframe:

# initiate plot
p <- ggplot(data=pix.df, aes(x=pix.df$time)) + ylab("VV")
# make sequence of pixel positions to plot
ind <- seq(2,9263, 67)
# loop
for (i in ind) { p <- p + geom_line(aes_string(y = pix.df[,i]), alpha = 0.2) }
p

or the first X columns, meaning a profile through the image:

p <- ggplot(data=pix.df, aes(x=pix.df$time)) + ylab("VV")
ind <- seq(2,60, 1)
for (i in ind) { p <- p + geom_line(aes_string(y = pix.df[,i]), alpha=0.3) }
p

Using Ratio to Estimate How Similar VV and VH Are

Instead of making a dataframe containing VV backscatter for each timestep, a VV/VH ratio is calculated per pixel, per timestep.

VV and VH seem to grow more similar during the course of the year. Surprising unison in first less and then more similarity is exhibited around the start of april.

Plotting Standard Deviation of VV and VH

Sentinel 2 data

NDVI_0204 <- read_stars("./S2_data/NDVI_02042018_clip.tif")
MDWI_3_8 <- read_stars("./S2_data/MDWI_3_8_02042018_clip.tif")
MDWI_8_12 <- read_stars("./S2_data/MDWI_8_12_02042018_clip.tif")
shape_NDVI <- st_transform(shape, crs=st_crs(NDVI_0204))
NDVI_0204_clip <- st_as_stars(NDVI_0204[shape_NDVI[124,]])
MDWI_3_8 <- st_as_stars(MDWI_3_8[shape_NDVI[124,]])
MDWI_8_12 <- st_as_stars(MDWI_8_12[shape_NDVI[124,]])
# plot(comb_split[1,,,9], main="VV 09.04.")
# plot(NDVI_0204_clip, main="NDVI 02.04.")
# plot(comb_split[2,,,9], main="VH 09.04.")

Example 1: April

plot(comb_split[1,,,9], main="VV 09.04.")
plot(comb_split[2,,,9], main="VH 09.04.")

Example 2: February

# NDVI from 01.02.2018, S1 from 27.01.
NDVI_0102 <- read_stars("./S2_data/NDVI_01022018_clip.tif")
NDVI_0102_clip <- st_as_stars(NDVI_0102[shape_NDVI[124,]])
plot(NDVI_0102_clip, main="NDVI 01.02.")
plot(comb_split[1,,,3], main="VV 27.01.")
plot(comb_split[2,,,3], main="VH 27.01.")

Time Series

Example 1: July

dates[16:18]
## [1] "2018-07-02" "2018-07-14" "2018-07-26"

Normalized Change Ratio

Normalized Difference VV and VH

Sentinel Playground screenshot (NDWI - plant moisture) in clarification of the dark spot in the middle left part of the polygon. NDVI shows slightly less Vegetation. Visual shows a slightly brighter patch. -> not water. sentinel playground screenshot - 16.07. / MDWI from bands 8A and 11

Example 2: February

Normalized Change Ratio